#include "navierstokes.hpp"

// variational formulation language
#include <life/lifevf/vf.hpp>

namespace Life
{


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J) //, vector_ptrtype& R)
{

    using namespace Life::vf;

    mesh_ptrtype mesh = Xh->mesh();

    element_fluid_type U( Xh, "u" ); U = *X;
    element_fluid_type V( Xh, "v" );
    element_fluid_0_type u = U.template element<0>();
    element_fluid_0_type v = V.template element<0>();
    element_fluid_1_type p = U.template element<1>();
    element_fluid_1_type q = V.template element<1>();

    //-----------
    double EE = 21*1e5;
    double sigma = 0.28;
    double mu = EE/(2*(1+sigma));
    double lambda = EE*sigma/((1+sigma)*(1-2*sigma));

    AUTO( Iddd , oneX()*trans(oneX()) + oneY()*trans(oneY()) );
    AUTO( Fv , Iddd + gradv(/*depl*/*M_aux) );
    AUTO( Ev, 0.5*(gradv(*M_aux)+trans(gradv(*M_aux)) ) + 0.5*trans(gradv(*M_aux))*gradv(*M_aux) );
    AUTO( Sv, lambda*trace(Ev)*Iddd + 2*mu*Ev );
    //-----------

   //--------------------------------------------------------------------------------------------------//

    AUTO( deft, 0.5*(gradt(u)+trans(gradt(u)) ) );
    AUTO( def, 0.5*(grad(v)+trans(grad(v)) ) );

    //Matrice Identitée
    AUTO (Id , oneX()*trans(oneX()) + oneY()*trans(oneY()) );

    // Tenseur des contraintes (trial)
    AUTO(Sigmat, -idt(p)*Id + 2*mu*deft );

    // Tenseur des contraintes (test)
    AUTO(Sigma, -id(q)*Id + 2*mu*def );

    //force volumique :
    //AUTO(f,vec(cst(0.),cst(0.)/*-rho*cst(9.81)*/ ) );
    AUTO(f,-rho*cst(9.81)*vec(cst(0.),cst(1.)) );

    // Vitesse d'entree du fluide
   // AUTO(Vit_Entree,Vitesse_Entree*( Py()*(1.0-Py()) )*oneX() );

    //--------------------------------------------------------------------------------------------------//
    //--------------------------------------------------------------------------------------------------//
    //--------------------------------------------------------------------------------------------------//
    //terme non lineaire
    form2( Xh,Xh, J , _init=true )  =
        integrate (elements(mesh), _Q<3*Order>(),
                   rho*trans(id(v))*(gradv(u))*idt(u));
    form2( Xh,Xh, J )  +=
        integrate (elements(mesh), _Q<3*Order>(),
                   rho*trans(id(v))*(gradt(u)*idv(u)) );

    form2( Xh,Xh, J )  +=
        integrate (elements(mesh), _Q<3*Order>(),
                   -rho*trans(idv(*M_mesh_velocity))*(gradt(u))*id(v) );


     // sigma : grad(v) sur Omega
    form2( Xh, Xh, J ) +=
        integrate( elements(mesh), _Q<2*(Order)>(),
            + trace(Sigmat*trans(grad(v))) );

    // terme assurant l'incompressibilité
    form2( Xh, Xh, J ) +=
        integrate( elements(mesh), _Q<2*(Order-1)>(),
            + divt(u)*id(q) );

    //--------------------------------------------------------------------------------------------------//
#if 1
    // terme de bord classique
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                   - trans(Sigmat*N())*id(v) );
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                   - trans(Sigmat*N())*id(v) );
    /*
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                   - trans(Sigmat*N())*id(v) );
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                   - trans(Sigmat*N())*id(v) );
    */
    //terme de bord weak boundary
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                   - trans(Sigma*N())*idt(u)
                   + penalbc*trans(idt(u))*id(v)/hFace() );
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                   - trans(Sigma*N())*idt(u)
                   + penalbc*trans(idt(u))*id(v)/hFace() );
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                   //- trans(Sigma*N())*idt(u)
                   + penalbc*trans(idt(u))*id(v)/hFace() );
    form2( Xh, Xh, J ) +=
        integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                   //- trans(Sigma*N())*idt(u)
                   + penalbc*trans(idt(u))*id(v)/hFace() );

#endif

    form2( Xh, Xh, J ) +=
        integrate( elements(mesh), _Q< 2*Order >(),
            + trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );
    //}

#if 0
 int stabilisation = this->application()->vm()["stab"].template as<int>();
if (stabilisation==1) {
   //terme de stabilisation
    form2( Xh, Xh, J ) +=
        integrate( internalfaces(mesh), _Q< 3*Order >(),
            + 0.025*vf::pow(hFace(),2.0)*( trans(idt(u))*N() )*(trans(jump(grad(u)))*jump(grad(v)) )  );
            //+ 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
//            + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
//( (trans(idt(u))*N())*( trans(id(u))*N() )) );

 form2( Xh, Xh, J ) +=
        integrate( elements(mesh), _Q< 3*Order >(),
 + 0.000*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
}
#endif


     J->close();

}

} // Life


// instantiation
template class Life::NavierStokes<2,2,Life::Simplex>;
//template class Life::NavierStokes<3,2,Life::Simplex>;

